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different integration methods. The numerical schemes we considered include spectral methods with 
different strategies for dealiasing and two variants of finite difference methods. Based on this com- 
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on the necessary numerical resolution are given. This estimate is based on analyzing the scaling 
behavior similar to the procedure in critical phenomena and present simulations are put into per- 
spective. 
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I. INTRODUCTION 

The question, whether the incompressible Euler equa- 
tions develop singularities in finite time starting from 
smooth initial conditions, remains an outstanding open 
problem in applied mathematics. Although substantial 
progress has been made in recent years using a more ge- 
ometrical viewpoint [H HI |3l [H [5] , it is yet not clear from 
numerical simulations, whether the assumptions of the 
theorems for non-blow up are fulfilled for flows evolving 
from simple smooth initial conditions. Singular struc- 
tures, evolving in finite time or simply "fast enough", 
may play a similar role as shock-like structures in com- 
pressible flows, providing structures which dominate the 
energy dissipation even in the non- viscous situation (see 
Eyink [6l [71 [8] and references therein) . 

In this paper, we study a Kida-Pelz like flow with dif- 
ferent numerical schemes: spectral methods with differ- 
ent strategies of dealiasing (this extends the study of Hou 
and Li [9 and confirms their results), two finite differ- 
ence methods and a finite volume method. Studying 
the structures of vorticity, it turns out that the differ- 
ences between the various methods of dealiasing are more 
pronounced than between the spectral methods and the 
finite difference/ volume methods. This result suggests 
that resolving the vorticity structures is more important 
than the order of the numerical scheme. It also justifies 
the use of finite difference/ volume methods in adaptive 
mesh refinement (AMR) simulations to resolve the vor- 
ticity structures. 

Using AMR simulations up to an effective resolution 
of 4096^ mesh points and comparing the results to lower 
resolution runs, we observe that the standard way of pre- 
senting a l/|u?| plot in time may lead to misleading con- 
clusions. However, looking at normalized plots reveals 
the issue of numerical resolution in a convincing manner. 
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II. NUMERICAL SCHEMES 

In this section we compare spectral methods with dif- 
ferent dealiasing and finite difference/ volume methods. 



A. Spectral methods and dealiasing 

We use a standard spectral method where the time 
stepping is performed with a strongly stable 3rd-order 
Runge-Kutta method [10] in Fourier space and where 
nonlinearities are calculated in real-space. On Linux- 
clusters, the FFTW-library is used whereas the library 
P3DFFT [11] from the San Diego Supercomputer Center 
is used on the IBM Regatta series and on BlueGene/L. 

We use three ways of dealiasing the spectral data: 

1. Spherical mode truncation: this is used in turbu- 
lence simulations (Biskamp and Miiller [12 ). The 
spherical mode truncation puts a sphere of radius 
y in Fourier space and nullifies all modes outside 
this sphere. 

2. Standard 2/3 rule: same as above, but using a ra- 
dius of |y = y [13 . This is the most common 
way of dealiasing spectral data. 

3. High-order exponential cut-off: this method was 
introduced by Hou and Li [9 and consists of in- 
troducing a high-order exponential filter function 
p{k) = ex.-p{—a{\k\/N)'^) with a = 36 and m = 36. 



B. Finite difference/ volumes methods 

All presented finite difference/ volume methods are sec- 
ond order and use the same strongly stable 3rd-order 
Runge-Kutta method \W as used in the spectral simula- 
tions. 

We implemented three different versions of real-space 
methods: 



1. Staggered grid formulation of Harlow and Welsh 
^14 : Normal components of the velocity are lo- 
cated at their respective cell faces and the pressure 
is defined at the cell centers. This allows a exact 
Hodge-decomposition such that no pressure oscil- 
lations occur. In addition, it conserves momentum 
and energy and could thus also been seen as a finite 
volume method. 

2. Vorticity formulation for AMR: From our previous 
AMR studies [151 IS] we know that the coarse-fine 
grid interpolations are very sensitive in the 3D Eu- 
ler simulations. As in the former simulations we 
choose to perform all data exchange and interpo- 
lation using the vorticity = V x u. Here, the 
vorticity is defined at cell centers and we applied a 
tri-cubic interpolation for coarse-fine grid interpo- 
lation. Then, three Poisson equations are solved for 
the cell-centered vector Potential A and staggered 
values for the velocity u = V x A are obtained. 

3. Finite volume method: this method is similar to 
the former but a finite volume method [15l [TTl [18] 
is used instead of finite differences. 



C. Comparison 

We first compare the growth of the maximum vortic- 
ity according to the Beale-Kato-Majda result [19l|20] for 
all six numerical methods described above. The initial 
condition was chosen similar to Kida-Pelz 12 vortices 
[2TJ [22l [23] with a Gaussian shape for the vorticity dis- 
tribution. Resolution of all the spectral simulations were 
512^ mesh points (corresponding to the full domain) and 
in addition the Hou-Li exponential filtering was repeated 
with 1024^ mesh points. The finite difference/ volume 
simulations were performed with 512^ and 1024^ mesh 
points. The growth of max|u?| is shown in Fig. [l] All 
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FIG. 1: Growth of max|u7| for all implemented numerical 
schemes. 
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FIG. 2: Isosurface plots of vorticity. From top to bottom: 
spherical truncation, 2/3 rule, exponential filtering, Harlow- 
Welsh, vorticity formulation, 512*^ mesh points 
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FIG. 3: Energy spectra at time t = 0.5 for spectral and finite 
difference methods: a) spherical model truncation (512^), b) 
high-order exponential cut-off (512^), c) 2/3 rule (512^), d) 
high-order exponential cut-off (1024^), e) vorticity formula- 
tion (1024^), f) staggered grid formulation (1024^) 



simulation agree very well up to the time when the flow is 
underresolved. This is about t = 0.4 for the simulations 
using 512^ mesh points and t = 0.47 for the 1024^ runs. 
There is no particular criterion which simulation per- 
forms better once the simulation is underresolved. The 
very simple message from this comparison is: you just 
have to resolve the flow and this is more important than 
the order of the scheme. 

In order to display the differences and similarities of 
the various numerical methods, we used a "low resolu- 
tion" simulation with 512^ mesh points at a late time 
t = 0.5 where the flow is already underresolved. There- 
fore, we looked at very low levels (5% of the maximum 
vorticity) as suggested and done by Kerr and Hou 
and Li [25]. Due to the high symmetry of the flow, only 
1/8 of the total configuration is shown. (To get a better 
impression for the geometry of the vortices, see Fig. |4j 
which shows an isosurface of 70% of the peak vorticity.) 

The spherical truncation produces highly visible arti- 
facts due to heavy oscillations which grow to substantial 
values. This is mostly suppressed in the simulation us- 
ing the classical 2/3 rule and nearly vanishes for the high 
order exponential smoothing. Thus our comparison con- 
firms the analysis of Hou and Li [25] . Remarkable is the 
strong similarity of the real-space methods to the spectral 
simulation with high order exponential smoothing. This 
is especially visible in Fig. [3] which shows the energy 
spectrum for spectral and finite difference/ volume meth- 
ods at time t = 0.5. In the spectral schemes, the spher- 
ical truncation and the 2/3 rule show a strong Gibbs 
phenomena which is absent in the exponential filtering 
and the finite difference/ volume schemes. The Har low- 
Welsh method is slightly more dissipative than the vor- 
ticity formulation. From the comparison with the spec- 
tral schemes using exponential filtering and 1024^ mesh 
points it is safe to say that the finite difference schemes 



with an approximately 1.3 times larger resolution in each 
spatial direction perform equally well as the spectral code 
with exponential filtering. Thus, our conclusions of this 
comparison is that the differences in the simulation re- 
sults caused by the choice of the dealiasing method are 
larger than the difference to and between the real-space 
methods. Our finding thus confirms the viewpoint of Or- 
landi and Carnevale [26 and justifies the use of finite 
difference/ volume methods as integration scheme in an 
adaptive mesh refinement treatment. 




FIG. 4: Isosurface plot of max|u7| at 70% of maximum vor- 
ticity. Shown is also the trajectory of a particle moving to the 
position of maximum vorticity. 
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FIG. 5: Growth of vorticity along the Lagrangian trajectory 
(red) which ends near the point of maximum vorticity and a 
fitted exponential (green). 
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D. Lagrangian trajectories 

As pointed out in [3l H] , the Lagrangian treatment of 
vorticity amplification is closely related to the local geo- 
metric properties - like curvature and torsion - of vortex 
lines. In Fig. |4]the trajectory of a Lagrangian tracer par- 
ticle is shown. To obtain this trajectory, we first iden- 
tified the spatial position of the maximum vorticity at 
a late time of the simulation and then traced back the 
actual trajectory. Fig. |5] shows the temporal evolution 
of vorticity following this trajectory. A tendency to an 
exponential growth of vorticity along the trajectory is 
obvious. 



III. ADAPTIVE MESH REFINEMENT 
SIMULATIONS 

A. The framework racoon 

For the adaptive mesh refine calculations, we use our 
framework racoon [27] which is designed for massive par- 
allel computations and scales for hyperbolic systems lin- 
early up to 16384 processors on BlueGene BG/L. How- 
ever, for the incompressible Euler equations, the pres- 
sure resp. vector potential are solved using an adaptive 
multigrid method [28 ^ 29 which presently scales only up 
to 64 processors. Therefore, the present simulations are 
limited to an effective resolution of 4096^ mesh points. 
Parallelization and load-balancing is performed using a 
space-filling Hilbert curve [27] . 

Using the framework racoon and the vorticity formula- 
tion, we solve the incompressible Euler equations with an 
effective resolution of 4096^ mesh points. Fig. |6] shows a 
volume rendering of vorticity at the latest time t = 0.5 
including the adaptive meshes. Memory consumption is 
quite moderate using less than 80 GBytes. 




FIG. 6: Volume rendering of vorticity at time t = 0.5. 
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FIG. 7: Temporal evolution of 1/ max |u?|. 



B. Analyzing the growth of vorticity 

Looking at Fig. [7| which shows the time evolution of 
l/max|u;| it is tempting to identify a finite time sin- 
gularity. However, a more appropriate presentation is 
obtained plotting max|u;| x (to — t) where to is the ex- 
pected singularity time. This quantity should converge 
to a horizontal line in this plot if a singularity occurs 
in finite time. The time to = 0.638 is chosen in a way 
that this scaling is observed in the late phase of the sim- 
ulation while the numerics is still resolved. This is 
shown in Figs. [8] and the zoom in the inlet of this figure. 
Especially the zoom of the late phase of the simulation 
demonstrates, how sensitive the growth of vorticity de- 
pends on the numerical resolution and that conclusions 
drawn from underresolved simulations must be handled 
with care. 
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FIG. 8: Scaling of the growth of vorticity. Red: 1024^ mesh 
points, Blue: 2048^ mesh points, Green: 4096^ mesh points. 
The inlet shows the late phase of the simulation and highlights 
the importance of numerical resolution. 
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IV. CONCLUSIONS AND OUTLOOK 

We demonstrated the extreme sensitivity of the growth 
of vorticity on the numerical resolution. In order to gain 
further insight into the mechanism of vorticity amplifi- 
cation, future simulations should include the following 
analysis and diagnostics: i) If a finite time singularity is 
expected, then the blow-up time of vorticity must occur 
at the same time when the spatial position of maximum 
vorticity and maximum strain come together, ii) The 
Lagrangian viewpoint should be analyzed according to 
Deng, Hou and Xu [4 and Gibbon [3]. iii) Simulations 
should use initial conditions including the Kida-Pelz fiow 
[2T] and Bob Kerr's orthogonal tubes [30 . However, the 
shape of the initial vortex tube should be chosen in such 



a way that vortex shedding will not pollute the vorticity 
growth. For orthogonal vortex tubes this was achieved by 
Orlandi and Carnevale [26^ starting with Lamb dipoles. 
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